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SUMMARY 

High Rayleigh number convection in a rectangular cavity with insulated 
horizontal surfaces and differentially heated vertical walls has been analyzed 
for an arbitrary aspect ratio smaller than or equal to unity. Unlike previous 
analytical studies, a systematic method of solution based on linearization 
technique and analytical iteration procedure has been developed to obtain 
approximate closed-form solutions for a wide range of aspect ratios. The pre- 
dicted velocity and temperature fields are shown to be in excellent agreement 
with available experimental and numerical data. 


INTRODUCTION 

It has long been recognized that the growth of crystals is affected by 
the fluid dynamics and heat transport characteristics of the process. The flow 
caused by natural convection has profound influence on the interface shape, 
defect density, and crystal integrity. As pointed out by Ostrach (ref. 1), 
high Rayleigh number natural convection flows could still exist under typical 
conditions in a reduced gravity environment, thereby influencing the growth of 
crystals. It is, therefore, important to understand the detailed mechanism of 
convection and flow behavior in order to produce high quality crystals in both 
ground- and space-based systems. 

A configuration containing all the essential physics describing the gen- 
eral flow behavior and transport processes related to crystal growth techniques 
is that of a rectangular enclosure heated from the side (ref. 2). It is impor- 
tant to note that a central problem inherent to all confined convection situa- 
tions is that the outer (core) flow region can not be predicted a priori and 
depends on the boundary layer which, in turn, is influenced by the core 
(ref. 3). This coupling of the boundary layer and the core region constitutes 
the main source of difficulty in studying the flow behavior. This matter is 
not merely a subtlety for analysis but has equal significance for numerical 
and experimental studies. Because of the complex nature of the problem, numer- 
ical and experimental studies must be performed with a high degree of accu- 
racy. Often a prior knowledge of flow behavior such as the extent of boundary 
layers or location of secondary cells can improve the numerical and experimen- 
tal studies by allowing the choice of proper mesh sizes and finer measurement 
techniques. It is therefore apparent that an analytical treatment of this 
problem, in particular for high Rayleigh number cases where the boundary layers 
develop, is of considerable importance. Such a treatment (especially if a 
closed-form solution can be obtained) could greatly enhance the understanding 
and predict the underlying physics of the problem. However, as stated by 
Ostrach (ref. 3), most available analytical studies contain unrealistic assump- 
tions and estimates regarding the core configuration which result in solutions 



with limited or undefined validity. Moreover, because of special assumptions 
used in analyzing each particular problem, there seems to be a lack of cohe- 
siveness and systematic approach in treating this type of problem. 

Tichy (ref. 4) presents a brief summary of recent work done on high 
Rayleigh number convection in low aspect ratio cavities. He, himself, obtained 
analytical-numerical solutions for the low aspect ratio enclosures. Unlike 
most of the previous analytical studies which concentrate on the prediction of 
an average Nusselt number, Tichy's work deals with the considerably more com- 
plex task of determining the detailed behavior of the flow and temperature 
fields. He assumes nearly parallel core flow with an aspect ratio c << 1. 

His predicted velocities agree reasonably well near centerline regions. In 
regions away from the centerline, the velocities become less accurate and the 
no-slip condition on vertical walls is no longer satisfied. Results of other 
recent investigations include velocity field prediction for e = 1 by Reddy 
et al. (ref. 5) for which they use a numerical penalty finite element method 
of calculation, and the recent work by Kamotani et al. (ref. 6) in which exper- 
imental data for Nusselt numbers are obtained for a wide range of important 
parameters . 

To the best of the author's knowledge there has not yet been an analytical 
investigation for high Rayleigh number natural convection flows for an arbi- 
trary aspect ratio (0 < e < 1). Most of the reported work has been done on 
very low aspect ratio enclosures and only a few of these studies are concerned 
with prediction of velocity fields. Moreover, as mentioned earlier, these con- 
tain many ad hoc assumptions concerning the nature of the core flow, thereby 
limiting the validity of solutions. It is the purpose of this study to develop 
an approximate but systematic analytical method of solution predicting the 
velocity and temperature distributions in a rectangular cavity with arbitrary 
aspect ratio 0 < e < 1 for the boundary layer region where the Rayleigh 
number Ra >> 1 . 

In the following analysis a combination of a linearizing method and an 
analytical iteration procedure is used to obtain approximate closed-form solu- 
tions to the problem. Analytical iteration is a useful technique which does 
not appear to have been fully utilized in fluid mechanics problems such as the 
present one (ref. 7). When iteration is used in conduction with another 
approximate procedure such as a linearization performed on the governing equa- 
tions, it leads to a powerful but straightforward technique for obtaining 
approximate, often closed-form, solutions to the problem. This report shows 
that the solutions obtained for this rather complex problem are in excellent 
agreement with existing experimental and numerical data. 
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b b = Pr -1 

g acceleration due to gravity 
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Gr Grashof number, Gr = 

V 

H height of cavity 

L width of cavity 

-1 - 1 /? 

i a = -8.31c Gr 

M M(y) = 0. 1 37y + m 

Nu Nusselt number, (eq. (71)) 

Nu r Nusselt number referred to in other studies 
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a 

Ra Rayleigh number, Ra * Pr • Gr 
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temperature 

temperature of the cold wall 
temperature of the hot wall 
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dimensionless horizontal velocity 
dimensionless vertical velocity 
horizontal coordinate 
dimensionless horizontal coordinate, 
vertical coordinate 
dimensionless vertical coordinate, y 
thermal diffusivity 


X 

L 


Y 

H 


3 



3 

A 

S 
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boundary-layer thickness for a vertical plate 
boundary-layer thickness 
aspect ratio, e = |- 

T - T 

dimensionless temperature, 9 = — ^ — 
x 1 .2s(e- 1/s -e-' /2s ) 


M = 


b 3/ V 1/4 

8m 3 


v kinematic viscosity 

? stream function, (eq. (3)) 


u/ dimensionless stream function, u/ = — 


reference stream function, ^ = vGr 


- 1/2 


Superscripts: 

* stretched coordinates 

B boundary layer 

c core 

CB core boundary layer 

CC core core 

1 first order approximate solutions 

2 second order approximate solutions 
Subscripts: 

t thermal boundary layer thickness 

l uni t order variables 
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MATHEMATICAL FORMULATION 


Figure 1 shows a schematic diagram of the physical system. The length of 
the rectangular cavity is assumed to be much larger than the height and the 
width so that the problem is two-dimensional. The top and bottom surfaces are 
thermally insulated and the left and right walls are maintained at isothermal 
hot and cold temperatures respectively. The flow field is assumed to be steady, 
laminar, and quasi-incompressible. By quasi-incompressibi 1 ty it is meant that 
the variation in density due to temperature differences are neglected except 
where they modify the body force term in the momentum equation. This is gener- 
ally referred to as Boussinesq approximation and has been discussed in a formal 
way in reference 8. 

The stream-function form of governing equations are 




c 4 A- + 2c 2 -4*^ ♦ ^1 


8x 


? ? 4 

8x 3y 9y J 
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ef R l xx yy J f xy V 
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where the stream function is defined by 


u = |^ 

8y 


V = - 


8<t/ 

8x 


(3) 


The pressure term has been eliminated by cross-differentiation of the momentum 
equations . 

The above equations are also nondimensional ized by use of the following 
variables: 


x 







(4) 


Notice that the reference stream function <j/ R is left unspecified and will be 
determined by a proper force balance as appropriate to the problem. For 
Gr >> 1 , Pr ~ 1 the correct force balance is obtained by equating the coeffi- 
cient of inertia and buoyancy terms in the vorticity equation (refs. 9 and 10). 
This yields to the proper characteristic stream function given by 


Y r = 


v a /g7 


(5) 


By substituting equation (5) into equations (1) and (2) the final form of the 
system of equations to be solved becomes 
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where 


and 
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b = 
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The boundary conditions are 
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(7) 


( 8 ) 


NORMALIZED BOUNDARY LAYER EQUATIONS 


The normalized boundary-layer equations are derived by properly balancing 
forces within the regions where they act. For Prandtl numbers larger than or 
equal to unity, the thermal boundary-layer thickness is smaller than, or at 
most equal to, the momentum boundary-layer thickness. In this case viscous 
and bouyant forces can then be assumed to be of the same order of magnitude 
within the thermal boundary-layer region (refs. 11 and 12). Stretching the 

B B 

coordinate normal to the vertical wall and expanding <|» and 0 in asymptotic 
series (ref. 9) gives 


* 

X 



(9) 


B B 

Y = E ] f 1 + . . . 


(10) 


B 

0 



+ . . . 


( 11 ) 


Substituting into general equations (6) and (7) and equating the coefficients 
of convective and conductive terms as well as viscous and buoyancy terms deter- 
mine the unknown parameters 6^. and E] 
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S t = ER a- ,/4 = eb ,/ V /2 
E, = Pr- 3/4 Gr-' /4 = b 3/4 a 1/2 


( 12 ) 

(13) 


and the boundary-layer equations given by 
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and 
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e i** = * e ] y + y ©! * 
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(15) 


The outer flow region or core equations are obtained by letting a 0 in 
equations (6) and (7). This gives 


C / 2 c C \ C / 2 c C \ C 
♦x c ^xxy + ^yyy J - ^y( e ^xxx + *Nyy ) = ®x 


(16) 


C C C C 

0/0 - u » 0 =0 

T x y T y x 


(17) 


The system of equations (6) and (7) has now been broken into the boundary- 
layer and core equations (14) to (17). An analytical iterative technique is 
used to find approximate solutions which satisfy equations (14) to (17). It 
is important to note that in an iterative process the accuracy of the solution 
depends on the number of iterations; therefore, if a limited number of itera- 
tions are desired, the choice of the initial profiles becomes very important. 
Consequently, if the initial profile is chosen so that it is most compatible 
with the particular assumptions and physics of the problem then the resulting 
solution will be satisfactory after only a few iterations. In the present 
problem the highest temperature changes occur within the thermal boundary- 
layer region near the vertical walls. Therefore the buoyant forces are very 
strong within that region, whereas in the core region the bouyant forces are 
weak because of the absence of sharp temperature gradients. Since the buoy- 
ancy is the main driving force here, it is important to choose an initial pro- 
file such that the change of buoyancy in the x-direction is accounted for. 
Accordingly we let 


1 

t = C ] (x)(16y 
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8y 2 + 1) 


(18) 


1 

0 = Dq(x) + (x) 



- 3y) 


(19) 


These profiles satisfy the horizontal wall conditions at y = + 1/2. Substi- 
tuting equations (18) and (19) into equations (14) to (17) yields the follow- 
ing equations: 
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C, ( 1 6y H - 8y + 1) 


3D } (4y^ - 1) 


B 1 fB 1 g < 
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and 

‘C I r c" ci r c 

C 1 ( 1 6y 4 - 8y 2 + 1) 16e 2 C 1 (4y 3 - y) + 3840^] - [16C 1 (4y 3 - y)J 

C i i * q • c 1 c 1 
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( 22 ) 


C’ 

C ] ( 1 6y 4 - 8y 2 + 1 ) 


3D 1 (4y - 1>J = p 6C ] (4y - y)J [D Q + D ] (4y J - 3y)J 


C 1 C’ 


(23) 


Note, that based on the discussion presented above, the functions C](x), 

D()(x), and Di ( x ) have been broken into boundary-layer and core contributions, 
thus accounting for varied buoyant action in the x-direction. Assuming the 
existence of an overlap region wherein both the core and the boundary-layer 
expansions remain valid, the following matching conditions must hold (see 
ref. 9): 
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IB 

Y* 

x-*» 
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= 0 x + -1/2 
1C 

= *x + -1/2 


IB 

t* 
x 

In order to solve the nonlinear and coupled boundary-layer and core equations 
(eqs. (20) to (23)), the following boundary-layer profiles, which satisfy the 
left wall conditions, are assumed: 


1C 

* ^x 
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(25) 


C, - Mb 3/4 Gr- |/4 


sin mx - cos mx 


; ) e”“] 
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* 

cos mxe 


mx 


(26) 


B 
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= L 


cos mxe 


An identical procedure can be used for the right wall. 


(27) 


The above functional forms are systematically predicted by utilizing the 
wall conditions, along with a linearized solution to the nonlinear equations. 
(See ref. 9 for a detailed description of "the modified Oseen linearization 
technique.") The unknown constants are determined by requiring the assumed 
profiles (18) and (19) to satisfy the integral form of governing equations in 
the y-direction and across the vertical boundary layer. The core region in 
the x-direction satisfies the core equations (16) and (17) exactly. The core 
and boundary-layer stream function and temperature are then 
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mx 
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4 (4y 1 * 3 
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(30) 


0 = -1 .78m 4 (4y 3 - 3y) 


(31) 


where m is an unknown constant to be determined later. The total stream 
function and temperature can be constructed by using the composition rule 
(ref. 13). We then obtai n 

1 .3/4 r -1/4 r / * *1 

t = oj 1 + (sin mx - cos mx) e mx (4y^ - 1)^ (32) 

8irr L v J 


1 - cos mxe mx ] (4y 3 - 3y) (33) 


6 = j e 


* 

mx 


- 1 . 78 nT 


To estimate m, the boundary-layer equation (14) is integrated 
* * 

from x = 0 to x = A/S at y = 0. The boundary- layer thickness A at 
y = 0 is taken to be the same as the boundary-layer thickness on a flat plate 
at that location, (ref. 14); this yields 

m = -0.56 (34) 
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Equations (28) to (33) are approximate solutions satisfying the boundary condi- 
tions on the walls and the integral form of governing equations. Although it 
is not shown here (for brevity), the above solutions predict the average 
Nusselt number and temperature field in good agreement with numerical and 
experimental data. Because of the averaged nature of the first approximation, 
however, the velocity fields are only in qualitative agreement with experimen- 
tal and numerical data. Therefore it seems reasonable that a refinement of 
solutions through an analytical iterative process should yield solutions in 
full agreement with existing studies. This is indeed the case as will be seen 
by comparing the iterated solutions to available data. (See section Results 
and Discussion.) 


By proceeding in the same manner as explained in reference 7, the approxi- 
mate results obtained for the core stream function and temperature (eqs. (30) 
and (31)) are utilized in the original vorticity equation to obtain a more 

1C 1C 

accurate core solution. The core variables y and 9 obtained previously 

2C 2C 

are a function of y only. Therefore, the new core variables and 0 

1C 1C 

can be expected to deviate only by a small margin from y and 9 . Hence 
the new variables must satisfy the condition 


2C 2C 


2C 2C 
99 99 

9x 9y 


(35) 


Using condition (35) in the original vorticity equation (6) results in the 
general core vorticity equation 


2C 2C 2C 
e ’•'yyyy + ^x ^ yyy 


2C 2C 
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= 0 
xyy x 


(36) 


The core temperature obtained in the first approximation is now substituted in 
the right-hand side of equation (36) to yield 
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2C 2C 
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e ^yyyy + ^yyy ^y ^xyy 


(37) 
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Assuming f = X(x)Y(y), we obtain 

x ' = d (38) 

and 

aY* 1 ’ 1 + ed(YY’ 1 1 - Y'y") = 0 (39) 


Note that equation (39) has boundary-layer characteristics (i.e., as a -> 0 
the highest derivative term vanishes); hence all boundary conditions on hori- 
zontal walls cannot be satisfied. Equation (39) can therefore be broken up 
into core and boundary-layer equations 
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dY 


CB 2CB 
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where y is the usual stretched coordinate 



for y § 0 


(41) 


(42) 


and superscripts "CB" and "CC" refer to "core boundary layer" and "core 
core" respectively. The parameter s is defined 


s 



(43) 


Omitting the details, it can be shown that equation (40) has an approximate 
solution which is 


y s A^y^ + A 2 y + A^ (44) 

The solution to the core boundary-layer equation (41) can be obtained by using 
an integral method, which is customary for finding approximate solutions to 
boundary-layer equations. A suitable profile that satisfies the horizontal 
wall conditions is given by 


*v 

1 - e~ y J (45) 

where C is a constant to be determined later. The profile equation (45) 
will later be forced to satisfy the integral form of equation (41). For the 
core stream functions we obtain 



> CB 

Y = C(ax + 1) 


= (Six + 1 ) ( A 1 y' 1 + A 2 y + A^) 

Equations (46) and (47) must match at an intermediate region (ref. 13) as 
fol lows 
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The composition rule for stream functions gives 

C CB CC CB 

= + (|i — ( X ,°°) 

Therefore the total core stream function is 


2C 


= C(Hx + 1) 




( 50 ) 


(51) 


where X, is a constant to be determined. It is now necessary to satisfy the 
1 CB 

core boundary-layer equation (41) by y on an integral basis. In addition, 
the "centro-symmetry" property (ref. 10) requires that 


2C 

<|> (0,0) = 0 (52) 


2C 

In the same manner as Collatz (ref. 7), the value of t at x = 0 and 

1C 

y * 0 will be set equal to the value of t at x * 0 and y = 0. The above 
constraints are sufficient to determine all constants in equation (51). These 
are 


where 


and 

where 


V 2s(e~' /s - r ,/2s ) 


S = 0.12 

3/4 -1/4 

C = An = % A 

8 nr 


A = 


(l - e" 1/2s ) 2 + X ] ( 1 /2s) 2 


(53) 


From relation (43), the parameter a is 


a = 


8.31 


eGr 


1/2 


(54) 


Therefore when Gr >> 
uni ty, and we obtain 


and 0 1 x | < 1/2, 2.x is negligible with respect to 


2C 
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*2 
+ X,y 


( 55 ) 
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A more accurate boundary-layer velocity can now be obtained by modifying 

IB 

the boundary-layer solution 14 / of equation (28) to satisfy the asymptotic 

2C 

matching condition with the new core stream function ^ of equation (55); 
that is, let 


2B 

Y - Ap 



COS 




(56) 


where the parameter M is an unknown function of y to be specified later. 
Note that the core-boundary-layer velocity matchings are satisfied; that is 


2B 


2C / 


t ( ra ,y> = t - 2>y 


(57) 


and 
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2B 2C/ 

^(®,y) = ff- j.y 
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(58) 


2B 

Also note that t satisfies the horizontal wall conditions. Assuming the 

2 B 

simplest linear profile for M(y) and requiring that is satisfied by the 

integral form of the boundary-layer energy equation (15) yields 


M(y) = 0. 137y + m 


(59) 


Again, using the composite rule for determining the total stream function, we 
obtai n 
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from which the velocities can be derived 
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r *n 
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sin (0. 137y + m)x - cos 

1(0. 1 37y + m)x 1 


.(0 


.137y + m)x j(^e y - e^ y + X^yj + 0.274Apj x s i n j^( 0 . 1 37y + m)xj 


(0. 137y + m)x 


1 - e~ y j + X,y 


(62) 


In a similar manner, a better approximation can be obtained for the core 

2C 

temperature. The new temperature 6 must satisfy the general core energy 
equation given by 
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Because of horizontal wall conditions imposed on © , it must be of the form 
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0 = g(y) + f(x) 

Substituting equation (64) into equation (63) yields 
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and hence 


f (x) = const. = 3 
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Using the horizontal wall condition and employing the central symmetry property 
2C 

0 (0,0) = 0 yield 
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Guided by the boundary-layer solution for 0, we let 


“ i .« 

0 = 2 e + 


1 - cos(mx)e 


mx 
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2 

where the core-boundary-layer matching conditions are satisfied for e. 

By the composition principle, the total temperature is constructed and 

reads 
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Defining the average dimensionless heat transfer coefficient at the wall by 

, 1/2 


Nu = 


0. 


y— 1/2 


x— 1/2 


dy 


(71) 


we then obtain 


Nu 


m_ 

25 


(1 + 3) 


(72) 


The reference Nusselt number used for comparison purposes can be shown to be 
related to Nu by 


and hence 


Nu r = eNu 


(73) 


N Uy . = ff-(l + 3) 


(74) 


To obtain an estimate for 3, we require that 0y and 0y be of the same 
order of magnitude in the 0 i y ^ 1/2 range (see ref. 7); that is, 


This gives 
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With determination of the parameter |3, the solution to the problem is 
complete. The approximate, closed-form solutions derived for the velocity and 
temperature fields now can be used conveniently to obtain detailed information 
about the flow behavior and heat transfer characteristics of the problem. In 
the next section the present solutions are carefully examined and comparisons 
to existing experimental and numerical data are made. 


RESULTS AND DISCUSSION 

In order to test the validity and accuracy of the solutions, a number of 
cases were calculated to enable a comparison with established experimental and 
numerical results of other investigations. 

The vertical and horizontal velocities as given by equations (61) and 
(62) have been plotted and compared to several existing numerical and experi- 
mental data for different values of the aspect ratio. 

Figure 2 shows a plot of the end vertical velocity profiles for 
Ra = 3x10^, Pr = 1, e = 0.1, and y = +0.179. The present solution in both 
the y = +0.179 and y = -0.179 regions .agrees very well with numerical data. 
The fact that Tichy's analytical solution does not satisfy the no-slip condi- 
tions at vertical walls except at y = 0 is partially responsible for a less 
quantitative agreement between his analytical and numerical solutions. Note 
that in the present study all wall conditions are satisfied. 

Figure 3 compares the present analytical solutions of mid-core velocities 
with Tichy's numerical solutions (ref. 4). The agreement is good for all 
velocities. It is especially important to note that the horizontal boundary- 
layer character of the flow for very large values of Rayleigh numbers is 
evident both in numerical and analytical solutions. Because of the centro- 
symmetric property, only the top part of velocities are drawn. The Rayleigh 
number ranges from 10? to 10^ with Pr = 1 and e = 0.1. 

In figure 4 the vertical velocity along the centerline of the enclosure 
is plotted and compared with numerical data of reference 5 which were obtained 
by using a penalty finite element method. There is good agreement through the 
entire region. In this case Ra = 10^, Pr = 1 , and the aspect ratio e is 
equal to unity (i.e., a square cavity). This, once again, emphasizes the wide 
range of applicability of the present solution. 

The average wall Nusselt number is plotted in figure 5. At very high 
Rayleigh numbers, Ra ~ 1 0 9 , the average heat transfer coefficient is practi- 
cally independent of the aspect ratio e. For moderately high Rayleigh numbers 
of Ra < 10^ the impact of the aspect ratio e becomes more obvious however, 
and results in a splitting of the Nusselt number curves which correspond to 
low aspect ratio enclosures and enclosures with aspect ratios near unity. The 
recent experiments by Kamotani , et al . (ref. 6) clearly support the current 
predictions . 

Figure 6 shows plots of the core temperature measured at three different 
x-locations and for the low aspect ratio of e = 0.0625 (ref. 6). The solid 
line is the core temperature distribution derived in the present analysis. 

There is good agreement between the analytical and experimental temperature 
fields. 
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CONCLUSIONS 


The problem of high Rayleigh number natural convection for arbitrary 
aspect ratio, 0 < e § 1 , has been solved analytically. A combination of a 
linearizing method and an analytical iteration technique was used. This combi- 
nation of techniques ultimately led to closed-form approximate solutions for 
velocity and temperature fields (eqs. (61), (62), and (70)). There are two 
major differences between this work and all of the known previous efforts; one 
is the wide range of applicability for the aspect ratio (0 < e S 1), and the 
other is the fact that the method of solution has been systematic, with no a 
priori assumptions made regarding, for example, the core configurations. This 
approach is especially important for even more complicated problems such as the 
study of multispecies natural convection where clear experimental and numeri- 
cal data are not available for use in theoretical analysis. 
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FIGURE 1. - PHYSICAL SYSTEM. 
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FIGURE 2. - END VERTICAL VELOCITY PROFILES OF ANALYTICAL 
AND NUMERICAL RESULTS; Ra = 3x10 8 , y = ±0.179, 

E = 0.1, Pr = 1. 
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